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1  Introduction 


This  is  the  fifth  article  of  a  series  [13, 14,  15,  16]  devoted  to  the  construction  and  study  of 
the  so-called  Runge-Kutta  Discontinuous  Galerkin  (RKDG)  method.  The  RKDG  method 
is  a  method  devised  to  numerically  solve  the  initial  boundary  value  problem  associated  with 
the  conservation  law 

-h  div  f (u)  =  0,  in  flx(0,T),  (1.1) 

where  Q  C  R!^  and  u  —  (ui,  which  is  assumed  to  be  hyperbolic,  that  is,  f(u) 

is  assumed  to  be  such  that  any  real  combination  of  the  Jacobians  J2i=i  has  m  real 
eigenvalues  and  a  complete  set  of  eigenvectors.  In  this  paper,  we  continue  our  work  in 
[13,  14,  15,  16]  and  extend  (and  improve)  the  RKDG  method  to  the  case  of  multidimensional 
systems.  To  place  this  paper  under  a  proper  perspective,  we  first  discuss  the  work  done  in 
this  series  of  papers  and  papers  by  other  authors  which  has  been  prompted  by  the  remarkable 
compactness  and  parallelizability  of  the  RKDG  method  and  by  its  ability  to  easily  handle 
boundary  conditions  and  complicated  geometry. 

The  original  discontinuous  Galerkin  method  was  introduced  by  Reed  and  Hill  [31],  and 
analyzed  by  LeSaint  and  Raviart  [26],  Johnson  and  Pitkaranta  [25],  Richter  [32],  and  by  Pe¬ 
terson  [29] .  All  these  were  for  the  linear  equations.  Our  work  was  concentrated  on  treating 
nonlinear  equations,  which  call  for  different  techniques.  The  first  (one-dimensional)  RKDG 
method  was  introduced  in  [13]  by  combining  the  piecewise-linear  discontinuous  finite  ele¬ 
ments  used  for  the  space  discretization  of  one-dimensional  conservation  laws  by  Chavent 
and  Cockburn  [11]  with  one  of  the  explicit,  TVD  time  discretizations  developed  by  Shu  [34], 
and  Shu  and  Osher  [35,  36].  The  resulting  scheme  was  shown  to  be  formally  uniformly 
second-order  accurate  (a  fact  confirmed  by  numerical  experiments)  and  was  proven  to  be 
total  variation  diminishing  in  the  means  (TVDM).  Later,  in  [14],  the  RKDG  schemes  were 
defined  using  a  general  framework  that  allowed  piecewise  polynomials  of  degree  k  E  N 
approximate  solutions.  These  fully  explicit  schemes  were  proven  to  be  TVBM  (total  vari- 
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ation  bounded  in  the  means)  and  were  shown  to  be  formally  uniformly  {k  +  l)-th  order 
accurate,  facts  that  were  both  verified  numerically.  The  extension  of  the  RKDG  schemes 
to  one-dimensional  systems  was  carried  out  in  [15]  and  the  multidimensional  case  for  the 
scalar  conservation  law  was  treated  in  [16],  where  it  was  proven  that  for  some  fairly  gen¬ 
eral  triangulations,  the  approximate  solution  given  by  the  RKDG  method  satisfies  a  local 
maximum  principle  independently  of  the  the  degree  k.  A  projection,  or  generalized  “slope 
limiting” ,  was  constructed  which  enforced  the  above  maximum  principle  without  destroying 
the  formal  accuracy  of  the  method.  Theoretical  indications  that  the  method  is  uniformly 
{k  +  l)-th  order-accurate  when  polynomials  of  degree  k  are  used  were  given  and  numerical 
validation  of  this  claim  was  presented  for  piecewise-linear  approximations  fc  =  1  in  uniform 
grids  made  of  triangles.  The  case  k  =  2  was  worked  out  by  Hou  [23].  An  extension  of  the 
RKDG  method  to  the  two-dimensional  Euler  equations  of  gas  dynamics  was  carried  out  in 
[17];  piecewise-linear  approximations  were  used.  In  this  paper,  we  complete  and  improve  the 
work  started  in  [17]. 

In  related  work,  Atkins  and  Shu  [1]  studied  an  alternative  quadrature-free  implementation 
of  the  RKDG  method.  Bey  and  Oden  [8]  used  the  RKDG  method  with  arbitrary  quadrilat¬ 
erals  and  piecewise-linear  approximate  solutions,  to  solve  2D  Euler  equations.  Jiang  and  Shu 
[24]  proved  a  cell  entropy  inequality  for  the  square  entropy  for  arbitrary  order  of  accuracy 
and  arbitrary  triangulations,  without  using  the  nonlinear  limiters,  for  the  semidiscrete  (con¬ 
tinuous  in  time)  case.  This  also  implied  the  stability  of  the  method  for  nonlinear  shocked 
cases.  Lowrie,  Roe  and  van  Leer  [27]  studied  the  discontinuous  Galerkin  method  in  space 
and  time;  see  also  the  related  studies  previously  made  by  Bar-Yoseph  [2]  and  Bar-Yoseph 
and  Elata  [3]. 

The  important  issue  of  the  parallelizability  of  the  RKDG  method  has  been  explored  by 
several  authors.  Biswas,  Devine,  and  Flaherty  [9]  have  shown  that  the  RKDG  method  (with 
a  new,  interesting  limiter)  has  a  “solution  parallel  efficiency”  of  99  %  in  the  NCUBE/2- 
a  reflection  of  the  fact  that  the  RKDG  method  uses  only  the  information  of  immediate 
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neighbors  to  march  in  time.  These  authors  have  also  constructed  h-  and  p-adaptive  versions 
of  the  RKDG  method  with  remarkable  results;  see  also  the  application  to  the  Euler  equation 
of  gas  dynamics  by  deCougny  et  at  [19].  The  important  issue  of  “dynamic  load  balancing,” 
essential  for  adaptive  methods,  has  been  addressed  by  Devine  et  al.  [21],  by  Ozturan  et  al. 
[28],  and  by  Devine  et  al  [20]. 

The  effect  of  the  quality  of  the  approximation  of  curved  boundaries  on  the  quality  of  the 

approximate  solution  has  been  explored  in  a  recent  paper  by  Bassi  and  Rebay  [4];  in  this 
paper,  we  only  consider  computational  domains  with  Lipschitz  boundaries. 

Extensions  of  the  method  to  the  compressible  Navier  Stokes  equations  and  general  con¬ 
vection  diffusion  equations  can  be  found  in  Bassi  and  Rebay  [5]  and  Cockburn  and  Shu  [18], 
respectively. 

We  are  now  ready  to  give  a  detailed  description  of  the  contents  of  this  paper.  In  Section  2, 
we  give  a  general  formulation  of  the  RKDG  method  for  multidimensional  systems,  including 
the  discussion  on  slope  limiters.  Section  3  contains  the  algorithm  and  implementation  details, 
including  the  numerical  fluxes,  quadrature  rules,  d^ees  of  freedom,  and  slope  hmiters  of 
the  RKDG  method  for  both  piecewise-linear  and  piecewise-quadratic  approximations  in  both 
triangular  and  rectangular  elements.  In  Section  4,  we  present  several  test  problems  for  the 
two-dimensional  Euler  equations  of  gas  dynamics  intended  to  illustrate  the  effect  of  the 
degree  k  and  the  effect  of  the  use  of  triangles  or  rectangles  on  the  accuracy  of  the  method. 
Concluding  remarks  are  given  in  Section  5. 

2  Algorithm  formulation 

To  define  the  RKDG  method,  we  proceed  as  in  [16]. 

2.1  Space  discretization 

First,  we  discretize  (1.1)  in  space  using  the  discontinuous  Galerkin  method.  For  each  time  t  E 
[0,  T],  the  approximate  solution  Uh{t)  is  sought  in  the  finite  element  space  of  discontinuous 
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functions 

i4  =  e  L-(fi)  -.VhWe  V{K),  yK  6  %},  (2.1) 

where  7^  is  a  triangulation  of  the  domain  il  and  V {K)  is  the  so-called  local  space.  In  this 
paper,  V{K)  is  taken  to  be  P’^,  the  collection  of  polynomials  of  degree  k,  for  k  =  l  and  2. 
To  determine  the  approximate  solution  Uh{t),  we  need  the  weak  formulation  of  (1.1): 


^  [  u{x,t)v{x)dx+  Y,  [  f{u{x,t)) -Ue^K  v{x)dr  -  [  f{u{x,t))  ■  gr&dv{x)dx  =  0, 

dtJK  eearc-'® 

for  any  smooth  function  v{x).  Here  denotes  the  outward  unit  normal  to  the  edge  e. 

We  replace  the  integrals  by  quadrature  rules  as  follows 

f{u{x,t))  ■ne^KVh{x:)dr  «  ■  ^e,K  v{Xel)\e\,  (2.2) 

1=1 

r  ^ 

f{u{x,t))  ■  gTa.dv{x)dx  «  (2.3) 

Jk 

Then,  the  flux  {{u{x,  t))  •  n^,*:  is  replaced  by  the  numerical  flux  he,K{x,  t),  the  exact  solution 
u  is  replaced  by  the  approximate  solution  Uh,  and  the  test  function  v  hyvneV  {K),  resulting 
in  the  following  scheme: 


Uh{t  =  0)  =  Pvi.  (no), 

d  r  ^ 

—  /  Uh{x,t)vh{x)dx  +  E  E  L0ihe,K{3:el,t)v{Xel)\e\ 

e€dK  1=1 

M 

-Y^j^i'^h{xKj,t))-graidvh{xKj)\K\=0,  yvheV{K),  ^/K  eTh.  (2.4) 
i=i 

The  operator  is,  for  example,  the  standard  L^-projection  into  the  finite  element  space 
Vh. 

The  value  of  the  numerical  flux  at  the  point  (x,  t),  he^K-(x,  f),  where  x  belongs  to  the  edge 
e  of  the  boundary  of  the  element  K,  depends  on  the  two  values  of  the  approximate  sohition 
at  {x,t).  One  is  the  value  obtained  from  the  mterior  of  the  element  K,  namely, 

f)  =  Yim.  Uh{y,t), 

y-^x,  yeK 
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and  the  other  is  the  value  obtained  from  the  exterior  of  the  element  K,  namely, 


7/1  (.«,*),  if  .X  6  (9Q, 

\imy^a:,y^KUh{y,t),  Otherwise. 


The  discrete  boundary  values,  7/1,  are  the  L^-projection  of  the  exact  boundary  data  7  into 
the  finite  element  space  obtained  by  taking  the  traces  of  the  elements  of  Vh  into  90. 

The  numerical  flux  is  deflned  as  he, /<:(•*)  t)  =  he,K{uh{x^'^^^^\  t),  t))  where  he^K 

is  any  two-pK>int  Lipschitz  flux,  which  is  monotone  in  the  scalar  case  and  is  an  exact  or 
approximate  Riemann  solver  in  the  system  case.  It  is  also  consistent  with  f('ti)  •  that 
is, 

he,K{u,u)  =  f{u)-  ne,K, 

and  conservative,  that  is. 


An  example  is  the  following  (local)  Lax-Priedrichs  flux 

he,K{o,,  b)  =  ^  [f(a)  •  rie./f  +  f(h)  •  rie^K  -  ae,K  {h- a)],  (2.5) 

where  ae,K  is  an  estimate  of  the  biggest  eigenvalue  of  the  Jacobian  ^f(w/i(x,  f))  ■  for 
(x,  t)  in  a  neighborhood  of  the  edge  e. 

It  is  convenient  to  take  the  local  spaces  V  {K)  to  be  the  space  of  polynomials  of  total 
degree  smaller  or  equal  to  k,  P^{K)\  in  this  case,  we  denote  I4  by  V^.  (Note  that  this 
choice  is  possible  regardless  of  the  shape  of  the  elements  K  since  the  functions  in  I4  are 
discontinuous.)  There  are  two  reasons  for  this  choice.  First,  if  the  local  space  V (K)  includes 
P'^{K),  it  is  possible  to  find  {fc+l)-th  order  accurate  approximations  in  V {K)  to  any  function 
in  Second,  if  V{K)  consists  of  polynomials  only  and  does  not  include  F*^+^(Rr), 

it  is  not  possible  to  find  {k  +  2)-th  order  accurate  approximations  in  V [K)  to  hmctions  in 
see  [12]. 

Moreover,  if  14  includes  ,  the  approximation  to  divf(n)  provided  by  the  above  space 
discretization  is  {k  +  l)-th  accurate  for  sufficiently  smooth  u,  provided  that  the  quadrature 
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rules  for  the  edges  of  the  elements,  (2.2),  are  exact  for  polynomials  of  degree  2fc  + 1,  and  the 
quadrature  rules  for  the  interior  of  the  elements,  (2.3),  are  exact  for  polynomials  of  degree 
2k,  see  [16],  Proposition  2.1.  It  is  thus  reasonable  to  expect  that  the  resulting  scheme  gives 
an  (k  +  l)-th  order  accurate  approximation  when  the  exact  solution  is  smooth  enough. 

For  the  choice  14  =  Vh  quadrature  rules  over  the  edges  exact  for  constants,  the 
resulting  scheme  is  nothing  but  a  finite  volume,  monotone  scheme  in  the  scalar  case.  Thus, 
the  discretization  by  the  Discontinuous  Galerkin  method  can  be  considered  as  a  high-order 
accurate  extension  of  finite  volume,  monotone  schemes. 

2.2  Time  discretization 

The  equations  defining  the  approximate  solution  can  be  rewritten  in  ODE  form  as  dt'^h  = 
Lh(uh,'ih)  after  inverting  the  “mass”  matrix.  Since  the  functions  of  14  are  discontinuous,  the 
“mass”  matrix  is  block-diagonal  and  the  blocks,  whose  orders  are  equal  to  the  dimensions 
of  the  local  spaces  V{K),  can  be  easily  inverted  by  hand.  If  a  locally  orthogonal  basis  is 
chosen,  the  mass  matrix  is  diagonal. 

If  we  are  using  a  finite  element  space  14  included  in  ,  we  would  like  to  discretize  in 
time  the  above  system  of  ODEs  with  a  method  that  is  at  least  {k  +  l)-th  order  accurate. 
To  do  that,  we  use  the  TVD  Riinge  Kutta  time  discretization  introduced  in  [34,  35].  Thus, 
if  is  a  partition  of  [0,T|  and  At"  =  -  r,n  =  0,  ...,N  -  1,  our  time-marching 

algorithm  reads  as  follows: 

•  Set  ul  =  Py^(uo); 

•  For  n  =  0, ...,  iV  —  1  compute  as  follows: 

1.  set  =  ul', 

2.  for  i  =  1, ...,  fc  -f- 1  compute  the  intermediate  functions: 

=  I  ^  +  ft;  Ar  L,(n«,  7^^”  +  ftAt”)) 

u=o 
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3. 


set 


=  u 


(fc+1) 

h 


Note  that  this  method  is  very  easy  to  code  since  only  a  subroutine  defining  Lh{uh,  7/i(*)) 
is  needed.  In  this  paper,  we  tise  the  second  order  and  third  order  accurate  Runge-Kutta  time 
discretizations  listed  below  in  Table  2.1,  for  piecewise  linear  and  piecewise  quadratic 
finite  element  approximations,  respectively. 


2.3  The  local  slope  limiting 

In  the  case  in  which  piecewise-constant  approximations  are  considered,  that  is,  when  14  = 
14°,  the  artificial  viscosity  that  the  numerical  flux  introduces  in  the  scheme,  due  to  upwind- 
ing,  is  enough  to  render  the  scheme  stable.  However,  when  the  local  spaces  are  richer,  the 
stabfiizing  influence  of  the  numerical  fluxes  is  not  enough  to  guarantee  the  absence  of  spu¬ 
rious  oscillations.  To  enhance  the  stability  of  the  method  and  eliminate  possible  spurious 
oscillations  in  the  approximate  solution,  a  local  slope  limiting  operator  Afl/i  is  introduced  in 
the  time-marching  algorithm  as  follows; 

•  Set  ul  =  Allft  Pvfc(«o); 

•  For  n  =  0, ...,  A/’  —  1  compute  as  follows; 

1.  set 
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2.  for  i  =  1, A:  +  1  compute  the  intermediate  functions; 

uf  =  An,  +  diAn)]  ; 

U=o  i 

3.  set 

Theoretical  studies  of  the  operator  AH,  can  be  found  in  [14]  for  the  one  dimensional  case, 
and  in  [16]  for  the  multidimensional  case.  Guided  by  these  results,  we  tise  in  this  paper  very 
simple,  practical,  and  effective  slope  limiting  operators  AH,.  To  compute  An,«,,  we  rely  on 
the  assumpUon  that  spurious  oscillations  are  present  in  Uh  only  if  they  are  present  in  its 
part  ul,  which  is  its  L^-projection  into  the  space  of  piecewise  linear  functions  V/^.  Thus,  if 
they  are  not  present  in  ul,  i.e.,  if 

ul  =  Anhul, 

then  we  assume  that  they  are  not  present  in  u,  and  hence  do  not  do  any  limiting: 


ATIhUh  =  Uh. 


On  the  other  hand,  if  spurious  oscillations  are  present  in  the  P^  part  of  the  solution  ul,  i.e., 
if 


U/^  An,  U,, 

then  we  chop  off  the  higher  order  part  of  the  numerical  solution,  and  limit  the  remaining  P^ 
part; 

An,u,  =  An,n^. 


In  this  way,  in  order  to  define  AH,  for  arbitrary  space  14,  we  only  need  to  actually  define  it 
for  piecewise  linear  functions  V^-  The  exact  way  to  do  that,  both  for  the  triangular  elements 
and  for  the  rectangular  elements,  will  be  discussed  in  the  next  section. 


3  Algorithm  and  implementation  details 

In  this  section  we  give  the  algorithm  and  implementation  details,  including  numerical  fluxes, 
quadrature  rules,  degrees  of  freedom,  fluxes,  and  limiters  of  the  RKDG  method  for  both 


piecewise-linear  and  piecewise- quadratic  approximations  in  both  triangular  and  rectangular 
elements. 

3.1  Fluxes 

For  the  numerical  flux  needed  in  (2.4),  we  use  the  simple  Lax-Eriedrichs  flux  (2.5): 

(a,  b)  =  ^  [f(a)  ■  rie^K  +  f(/^)  •  -  0!e,K  (b  -a)]- 

The  numerical  viscosity  constant  ae,K  should  be  an  estimate  of  the  biggest  eigenvalue  of  the 
Jacobian  ^f(uh(x,  t))  ■  rie^K  for  {x,  t)  in  a  neighborhood  of  the  edge  e.  For  the  triangular 
elements,  we  have  used  the  local  Lax-Friedrichs  recipe; 

•  Take  ae,K  to  be  the  larger  one  of  the  largest  eigenvahie  (in  absolute  value)  of 

ne,K  and  that  of  •  Ue^K,  where  uk  and  Uk>  are  the  means  of  the  numerical 

solution  in  the  elements  K  and  K'  sharing  the  edge  e. 

For  the  rectangular  elements,  we  have  used  both  the  local  Lax-Friedrichs  recipe  (in  Ex¬ 
amples  4.1  and  4.2)  and  the  global  Lax-Friedrichs  recipe  (in  Example  4.3); 

•  Take  o-e./c  to  be  the  largest  of  the  largest  eigenvalue  (in  absolute  value)  of  ■^f{uK")-ne^K, 
where  Uk"  is  the  mean  of  the  numerical  solution  in  the  element  K",  which  runs  over 
all  elements  on  the  same  line  (horizontally  or  vertically,  depending  on  the  direction  of 
77,6, k)  with  K  and  K'  sharing  the  edge  e. 

Usually,  the  global  Lax-Friedrichs  recipe  is  more  dissipative,  but  is  more  robust,  than 
the  local  Lax-Friedrichs  recipe,  especially  for  problems  involving  low  velocities  and  low  den¬ 
sity/pressure  near  wall  boundaries.  There  are  recipes  in  between  the  two,  such  as  taking  the 
maxiTiniim  over  several  neighboring  elements  in  obtaining  oee,^:,  but  we  have  not  iised  them 
in  this  paper. 
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3.2  Quadrature  rules 

According  to  the  analysis,  the  quadrature  rules  for  the  edges  of  the  elements,  (2.2),  must 
be  exact  for  polynomials  of  degree  2A:  +  1,  and  the  quadrature  rules  for  the  interior  of  the 
elements,  (2.3),  must  be  exact  for  polynomials  of  degree  2k,  if  methods  are  used.  Here  we 
discuss  the  quadrature  points  used  for  and  in  the  triangular  and  rectangular  element 
cases. 

3.2.1  The  rectangular  elements 

For  the  edge  integral,  we  use  the  following  two  point  Gaussian  rule 


for  the  F^  case,  and  the  following  three  point  Gaussian  rule 

£9W<ix«®[5(-3)+9(®)]+®9(0),  (3.2) 

for  the  F^  case,  suitably  scaled  to  the  relevant  intervals. 

For  the  interior  of  the  elements,  we  could  use  a  tensor  product  of  (3.1),  with  4  quadrature 
points,  for  the  F^  case.  But  to  save  cost,  we  “recycle”  the  values  of  the  fluxes  at  the  element 
boundaries,  and  only  add  one  new  quadrature  point  in  the  middle  of  the  element.  The 
quadrature  nile  is  thus: 

f_j\(x,y)dxdy  «  i[!,  ->.“75) +»  "TS'-O""® 

71'')  -^,iy  +29(0,0). 

For  the  F^  case,  we  use  a  tensor  product  of  (3.2),  with  9  quadrature  points. 

3.2.2  The  triangular  elements 

For  the  edge  integral,  we  use  the  same  two  point  or  three  point  Gaussian  quadratures  as  in 
the  rectangular  case,  (3.1)  and  (3.2),  for  the  F^  and  F^  cases,  respectively. 
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For  the  interior  integrals  (2.3),  we  use  the  three  mid-point  rule 

f  9{x,y)dxdy  fa  > 

•'!<  i=i 

where  m.j  are  the  mid-points  of  the  edges,  for  the  case.  For  the  case,  we  use  a  7  point 
quadrature  rule  which  is  exact  for  polynomials  of  degree  5  over  triangles,  given  in  Table  A.4, 
on  page  343  of  [10]. 

3.3  Basis  and  degrees  of  freedom 

We  emphasize  that  the  choice  of  basis  and  degrees  of  freedom  does  not  affect  the  algorithm, 
as  it  is  completely  determined  by  the  choice  of  function  space  V{h)  in  (2.1),  the  numerical 
fluxes  in  (2.4),  the  quadrature  rules,  the  slope  limiting,  and  the  time  discretization.  However, 
a  suitable  choice  of  basis  and  degrees  of  freedom  may  simplify  the  implementation  and 
calculation. 

3.3.1  The  rectangular  elements 

For  the  P^  case,  we  use  the  following  expression  for  the  approximate  solution  Uh{x,  y,  t) 
inside  the  rectangular  element  x  [yj_i.,yj^i]: 

Uh{x,  y,  t)  =  u{t)  +  u^{t)(l>i{x)  +  Uy{t)i}j{y)  (3.3) 


where 


and 


<f>i{x) 


X  —  Xj 
Axi/2' 


i’jiy)  = 


y  -  yj 

^yj/2^ 


Axi  =  Xi^i  -Xi_i  Ayj  =  yj^  -  y^.i- 


The  degrees  of  freedoms,  to  be  evolved  in  time,  are  then 


(3.4) 


'u(t),  Ux{t))  'Uj,(t)- 

Here  we  have  omitted  the  subscripts  ij  these  degrees  of  freedom  should  have,  to  indicate 
that  they  belong  to  the  element  ij  which  is  [x^_i ,  rj+i] 
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Notice  that  the  basis  functions 


1,  (l)i{x), 

are  orthogonal,  hence  the  local  mass  matrix  is  diagonal: 

M  =  AxiAyj  diag  (^1,  . 

For  the  case,  the  expression  for  the  approximate  solution  Uh{x,  y,  t)  inside  the  rectan¬ 
gular  element  [x^  X  [y^  i,yj+i]  is: 

Uh{x,  y,  t)  =  u{t)  +  u^{t)(l>i{x)  +  Uy{t)i)j{y)  -f  u^y{t)^i{x)'ipj{y) 

+u^,{t)  -  i)  +  K„(t)  (./.Jfo)  -  1)  .  (3.5) 

where  0i(.x)  and  'ipjiy)  are  defined  by  (3.4).  The  degrees  of  freedoms,  to  be  evolved  in  time, 
are 

Uyy{t). 

Again  the  basis  functions 

•1,  (l)i{x),  -ipj{y),  4>i{x)'^j{y), 

are  orthogonal,  hence  the  local  mass  matrix  is  diagonal: 

A  A  /  1  1  1  4  4  \ 

M  -  AxiAyj  diag  (^1,  ,  — j  . 

3.3.2  The  triangular  elements 

For  the  P^  case,  we  use  the  following  expression  for  the  approximate  solution  Uh{x,y,t) 
inside  the  triangle  K: 

3 

Uh{x,y,t)  ^'^Ui{t)(pi{x,y) 

i=l 

where  the  degrees  of  freedom  Ui{t)  are  values  of  the  numerical  solution  at  the  midpoints  of 
edges,  and  the  basis  function  (pi{x,y)  is  the  linear  function  which  takes  the  vahie  1  at  the 
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mid-point  of  the  f-th  edge,  and  the  vahie  0  at  the  mid-points  of  the  two  other  edges.  The 
mass  matrix  is  diagonal 

For  the  case,  we  use  the  following  expression  for  the  approximate  solution  Uh{x,y,t) 
inside  the  triangle  K: 

6 

Uh{x,y,t)  =  Y^Ui{t%{x,y) 
i=l 

where  the  degrees  of  freedom,  Ui{t),  are  vahies  of  the  numerical  sohition  at  the  three  mid¬ 
points  of  edges  and  the  three  vertices.  The  basis  function  ^i{x,y),  is  the  quadratic  function 
which  takes  the  value  1  at  the  point  i  of  the  six  points  mentioned  above  (the  three  midpoints 
of  edges  and  the  three  vertices),  and  the  value  0  at  the  remaining  five  points.  The  mass 
matrix  this  time  is  not  diagonal. 

3.4  Limiting 

We  construct  slope  limiting  operators  AUh  on  piecewise  linear  functions  Uh  in  such  a  way 
that  the  following  properties  are  satisfied; 

1.  Accuracy:  if  Uh  is  linear  then  kXihUh  =  u^. 

2.  Conservation  of  mass:  for  every  element  K  of  the  triangulation  7^,  we  have: 

/  ATihUh  =  /  Uh. 

Jk  Jk 

3.  Slope  limiting:  on  each  element  K  of  %,  the  gradient  of  AYihUh  is  not  bigger  than 
that  of  Uh- 

The  actual  form  of  the  slope  limiting  operators  is  closely  related  to  that  of  the  slope 
limiting  operators  studied  in  [14]  and  [16]. 
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3.4.1  The  rectangular  elements 


The  limiting  is  performed  on  and  Uy  in  (3.3),  using  the  differences  of  the  means.  For  a 
scalar  equation,  Ux  would  be  limited  (replaced)  by 

m  {ux,  Ui+i  j  -  Uij,  Uij  -  Uj-ij)  (3.6) 


ai,  if  jail  <  MAx'^- 

m(ai, ...,  am),  otherwise 


(3.7) 


s  mini  jail,  if  s  =  sign{ai) 
0,  otherwise 


sign{am)-, 


where  the  function  m  is  the  TVB  corrected  minmod  function  [33,  14]  defined  by 

ra  (fii, ...,  Om^ 

with  the  minmod  function  m  defined  by 
m(ai,  ...,am)  = 

The  TVB  correction  is  needed  to  avoid  unnecessary  limiting  near  smooth  extrema,  where  the 
quantity  Ux  or  Uy  is  on  the  order  of  0{Ax'^)  or  0{Ay^).  For  an  estimate  of  the  TVB  constant 
M  in  terms  of  the  second  derivatives  of  the  function,  see  [14].  Usually,  the  numerical  results 
are  not  sensitive  to  the  choice  of  iff  in  a  large  range.  In  all  the  calculations  in  this  paper  we 
take  M  to  be  50. 

Similarly,  Uy  is  limited  (replaced)  by 

m(uy,  Ui  j^i  aij,  Uij  '^i,j  i)- 


with  a  change  of  Ax  to  Ay  in  (3.7). 

For  systems,  we  perform  the  limiting  in  the  local  characteristic  variables.  To  limit  the 
vector  Ux  in  the  element  ij,  we  proceed  as  follows: 

•  Find  the  matrix  R  and  its  inverse  R~^,  which  diagonalize  the  Jacobian  evaluated  at 
the  mean  in  the  element  ij  in  the  a:-direction: 

ou 

where  A  is  a  diagonal  matrix  containing  the  eigenvalues  of  the  Jacobian.  Notice  that 
the  columns  of  R  are  the  right  eigenvectors  of  and  the  rows  of  are  the  left 

eigenvectors. 
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Figure  3.1:  Illustration  of  limiting. 

•  Transform  all  quantities  needed  for  limiting,  i.e.,  the  three  vectors  u^ij,  Ui+i,j  —  Wy  and 
Uij  —  Ui-i  j,  to  the  characteristic  fields.  This  is  achieved  by  left  multiplying  these  three 
vectors  by  i2“h 

•  Apply  the  scalar  limiter  (3.6)  to  each  of  the  components  of  the  transformed  vectors. 

•  The  result  is  transformed  back  to  the  original  space  by  left  multiplying  R  on  the  left. 

3.4.2  The  triangular  elements 

To  construct  the  slope  limiting  operators  for  triangular  elements,  we  proceed  as  follows.  We 
start  by  making  a  simple  observation.  Consider  the  triangles  in  Figure  3.1,  where  mi  is  the 
mid-point  of  the  edge  on  the  boundary  of  Ko  and  bi  denotes  the  barycenter  of  the  triangle 
Ki  for  i  =  0, 1,2,3. 

Since  we  have  that 

mi  -  bo  =  ai  {bi  -  bo)  -h  0:2  (^2  -  bo), 

for  some  nonnegative  coefficients  ai ,  0:2  which  depend  only  on  mi  and  the  geometry,  we  can 
write,  for  any  linear  function  Uh, 

Uhimi)  -  Uh{bo)  =  Q!i  {uh{bi)  -  Uh{bo))  +  0:2  {uh{b2)  -  Uh{bo)), 
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and  since 


we  have  that 


UKi  ~  I  I  [  '^h  —  Uh{bi), 
\Ki\  JKi 


i  =  0,l,2,3, 


Uh{mi,  Ko)  =  Uh{mi)  -  uk„  =  ai  (u/Ci  -  w/c,,)  +  “2  («K2  -  uk„)  =  Au{mi,Ko) 

Now,  we  are  ready  to  describe  the  slope  limiting.  Let  us  consider  a  piecewise  linear  function 
Uh,  and  let  ruj,  i  =  1, 2, 3  be  the  three  mid-points  of  the  edges  of  the  triangle  Kq.  We  then 
can  write,  for  {x,y)  6  Kq, 

3  3 

Uh{x,y)  =  J2'^h{mi)<Piix,y)  =  uk^  +  Y^Uh{'mi,  Kf))(pi{x,y). 

i=l  i=l 

To  compute  AUhUh,  we  first  compute  the  quantities 


Aj  =  m{uh(mi,  Ko),i'Au{mi,  Kq)), 

where  m  is  the  TVB  modified  minmod  function  defined  in  (3.7),  and  v  >  \.  We  take  i/  =  1.5 
in  oxir  mmierical  runs.  Then,  if  A*  =  0,  we  simply  set 

3 

AUhUh{x,  y)  =  uk„  +  X!  Ai  (pi{x,  y). 

i=l 

If  X)|=i  Aj  ^  0,  we  compute 

3  3 

po.s  =  ^max(0,  Aj),  ne^  =  ^max(0,  —  Aj), 

i=l  i=l 


and  set 


Then,  we  define 


where 


pos  )  neg 


AUhUh{x,  y)  =  uk„  + 


Aj  =  S'*' max(0,  Ai)  —  0  max(0,  —  Aj). 
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It  is  very  easy  to  see  that  this  slope  limiting  operator  satisfies  the  three  properties  listed 

above. 

For  systems,  we  perform  the  limiting  in  the  local  characteristic  variables.  To  limit  Aj,  we 
proceed  as  in  the  rectanguiar  case,  the  only  difference  being  that  we  work  with  the  following 
Jacobian 


rrii  -  bp 
\mi  -  bo 


4  Numerical  Results 


In  this  section  we  present  several  mimerical  results  obtained  with  the  and  (second 
and  third  order  accurate)  RKDG  methods  with  either  rectangular  or  triangular  elements. 
These  are  standard  test  problems  for  Euler  equations  of  compressible  gas  dynamics. 


Example  4.1.  Double  Mach  reflection  of  a  strong  shock.  This  problem  was  studied  exten¬ 
sively  in  Woodward  and  Colella  [37]  and  later  by  many  others.  We  use  exactly  the  same 
setup  as  in  [37],  namely,  a  Mach  10  shock  initially  making  a  60°  angle  with  a  reflecting  wall. 
The  undisturbed  air  ahead  of  the  shock  has  a  density  of  1.4  and  a  pressure  of  1. 

We  use  rectangular  elements  for  this  problem.  The  compxitational  domain  is  [0, 4]  x  [0, 1], 
as  in  [37].  The  reflecting  wall  lies  at  the  bottom  of  the  computational  domain  for  |  <  a;  <  4. 
Initially  a  right-moving  Mach  10  shock  is  positioned  at  x  =  =  0  and  makes  a  60°  angle 

with  the  x-axis.  For  the  bottom  boundary,  the  exact  post-shock  condition  is  imposed  for 
the  part  from  a;  =  0  to  .•r  =  |,  to  mimic  an  angled  wedge.  Reflective  boundary  condition  is 
used  for  the  rest.  At  the  top  boundary  of  our  computational  domain,  the  flow  values  are  set 
to  describe  the  exact  motion  of  the  Mach  10  shock.  Inflow/outflow  boundary  conditions  are 
used  for  the  left  and  right  boundaries.  The  resxilts  at  t  =  0.2  are  shown.  As  in  [37],  only  the 
results  in  [0, 3]  x  [0, 1]  are  displayed. 

Four  different  uniform  meshes  are  used:  240  x  60  elements  (Ax  =  Ay  =  ^);  480  x  120 
elements  (Ax  =  Ay  =  960  x  240  elements  (Ax  =  Ay  =  ^);  and  1920  x  480  elements 

(A,x  =  Ay  =  ^)-  The  density  is  plotted  in  Figure  4.1  for  the  P^  case  and  in  Figure  4.2 
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for  the  case.  In  all  the  plots,  we  use  30  contours  equally  distributed  from  p  =  1.3965  to 

p  =  22.682. 

It  is  not  easy  to  observe  any  significant  difference  between  the  P^  and  P^  results  in  these 
pictures.  However,  if  we  show  a  “blowed  up”  portion  around  the  double  Mach  region,  in 
Figure  4.3,  we  can  see  that  P^  with  Ax  =  Ay  =  has  qualitatively  the  same  resolution 
as  P^  with  Ax  =  Ay  =  -^  for  the  fine  details  of  the  complicated  structure  in  this  region. 
Notice  that  this  detailed  structure  is  of  physical  interest  and  was  studied  before  with  an 
adaptive  grid  calculation  in  [7].  P^  with  Ax  =  Ay  =  -^  gives  a  much  better  resolution  for 
these  structures  than  P^  with  the  same  number  of  elements. 

The  conclusion  here  is  that,  if  one  is  interested  in  such  fine  structures,  then  one  can  use 
the  third  order  scheme  P^  with  only  half  of  the  mesh  points  in  each  direction  as  in  Ph  This 
translates  to  a  reduction  of  a  factor  of  8  in  space-time  cells  for  2D  time  dependent  problems, 
and  will  more  than  off-set  the  increase  of  cost  per  cell  and  the  smaller  CFL  number  by  using 
the  higher  order  P^  method  (the  cpu  saving  for  this  problem  is  around  a  factor  of  2.1  in  our 
implementation).  This  saving  wiU  be  even  more  significant  for  3D. 

The  optimal  strategy,  of  course,  is  to  use  adaptivity  and  concentrate  cells  around  the 
interesting  region,  and/or  change  the  order  of  the  scheme  in  different  regions. 

Example  4.2.  Flow  past  a  forward  facing  step.  This  problem  was  again  studied  extensively 
in  Woodward  and  Colella  [37]  and  later  by  many  others.  The  setup  of  the  problem  is  the 
following;  a  right  going  Mach  3  uniform  flow  enters  a  wind  tunnel  of  1  unit  wide  and  3 
iinits  long.  The  step  is  0.2  units  high  and  is  located  0.6  units  from  the  left-hand  end  of 
the  tunnel.  The  problem  is  initialized  by  a  uniform,  right-going  Mach  3  flow.  Reflective 
boundary  conditions  are  applied  along  the  walls  of  the  tunnel  and  in-flow  and  out-flow 
boundary  conditions  are  applied  at  the  entrance  (left-hand  end)  and  the  exit  (right-hand 
end),  respectively.  The  results  at  t  =  4  are  shown. 

The  corner  of  the  step  is  a  singularity,  which  we  study  carefully  in  our  numerical  experi¬ 
ments.  Unlike  in  [37]  and  in  many  other  papers,  we  do  not  modify  our  scheme  near  the  corner 
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Figure  4,1:  Double  Mach  reflection  problem.  Second  order  results.  Density  p.  30  equally 
spaced  contour  lines  from  p  =  1.3965  to  p  =  22.682.  Mesh  refinement  study.  Prom  top  to 
bottom:  Aa:  Ay  =  and 
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Figure  4.2:  Double  Mach  reflection  problem.  Third  order  results.  Density  p.  30  eqiially 
spaced  contour  lines  from  p  =  1.3965  to  p  =  22.682.  Mesh  refinement  study.  Prom  top  to 
bottom:  Aa;  =  Ay  =  ^ ^ ,  and  . 
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0.5 


Rectangles  P2,  A  x  =  A  y  =  1/240 


Figure  4.3:  Double  Mach  reflection  problem.  Blowed-up  region  around  the  double  Mach 
stems.  Density  p.  Third  order  with  Ax  =  Ay  =  ^  (top);  second  order  P‘  with 
Ax  =  Ay  =  ^  (middle);  and  third  order  P^  with  A.x  =  Ay  =  (bottom). 
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in  any  way.  It  is  well  known  that  this  leads  to  an  erroneous  entropy  layer  at  the  downstream 
bottom  wall,  as  well  as  a  spurious  Mach  stem  at  the  bottom  wall.  However,  these  artifacts 
decrease  when  the  mesh  is  refined.  In  Figure  4.4,  second  order  results  using  rectangtilar 
elements  are  shown,  for  a  mesh  refinement  study  using  A.r  =  Ay  =  and  ^  as 

element  sizes.  We  can  clearly  see  the  improved  resolution  (especially  at  the  upper  slip  line 
from  the  triple  point)  and  decreased  artifacts  caused  by  the  corner,  with  decreased  element 
sizes.  In  Figure  4.5,  third  order  results  using  the  same  sequence  of  elements  are  shown. 
Comparing  with  the  results  in  Figure  4.4,  we  can  see  that  the  resolution  is  improved, 
especially  for  the  slip  line  issued  from  the  triple  point. 

In  order  to  verify  that  the  erroneoiis  entropy  layer  at  the  downstream  bottom  wall  and  the 
spurious  Mach  stem  at  the  bottom  wall  are  both  artifacts  caused  by  the  poor  resolution  of 
the  corner  singularity,  we  use  our  triangle  code  to  locally  refine  near  the  comer  progressively. 
A  sequence  of  such  triangulation  is  shown  in  Figure  4.6,  where  cr  is  the  ratio  between  the 
typical  size  of  the  triangles  near  the  corner  and  that  elsewhere.  The  resolution  of  the  meshes 
away  from  the  corner  is  roughly  equal  to  a  rectangular  element  case  of  A.x  =  Ay  =  ^,  i.e., 
the  top  pictures  in  Figures  4.4  and  4.5.  The  density  results  using  P^  and  these  triangiilations 
are  shown  in  Figtire  4.7,  those  using  are  shown  in  Figure  4.8.  We  can  see  that,  with  more 
triangles  concentrated  near  the  corner,  the  artifacts  gradually  decrease.  Notice  that  there  is 
a  strong  spurious  entropy  production  near  the  corner,  which  pollutes  the  flow  downstream. 
This  is  more  apparent  when  the  entropy  is  plotted  (pictures  not  shown  to  save  space).  With 
progressive  refinement  near  the  comer,  this  spurious  entropy  production  decreases. 

These  are  the  only  triangular  element  runs  we  present  in  this  paper.  We  can  see  that 
the  triangular  elements  can  give  results  of  the  same  resohition  quality  as  the  rectangular 
case,  with  roughly  the  same  mesh  density,  for  both  P^  and  P^.  We  do  observe,  however, 
that  a  positivity  correction  procedure  is  needed  for  the  triangular  element  runs  for  this  case. 
During  the  projection  of  the  linear  part,  we  check  whether  the  density  and  the  total  energy 
are  negative  at  the  three  mid-points  of  the  edges  of  K.  If  they  are,  further  limiting  is 
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performed  to  bring  them  to  in  a  conservative  way. 

Example  4.3.  Shock  passing  a  backward  facing  comer  (diffraction).  This  example  has  been 
used  in  [22]  and  [30],  see  also  the  experimental  results  in  [6],  The  setup  of  the  problem  is  the 
following:  the  computational  domain  is  the  union  of  [0, 1]  x  [6, 11]  and  [1, 13]  x  [0, 11];  the 
initial  condition  is  a  pure  right  moving  shock  of  Mach  =  5.09,  initially  located  at  a;  =  0.5 
and  6  <  y  <  11,  moving  into  an  undisturbed  air  ahead  of  the  shock  with  a  density  of  1.4 
and  a  pressure  of  1.  The  boundary  conditions  are  inflow  at  .a:  =  0, 6  <  y  <  11,  outflow  at 
a;  =  13, 0  <  y  <  11,  reflective  at0<a;<  l,y  =  6  and  at  .x  =  1, 0  <  y  <  6,  and  Neumann 
at  1  <  a;  <  13,  y  =  0  and  at  0  <  .x  <  13,  y  =  11.  No  special  treatment  is  done  at  the  corner 
which  is  a  singularity  of  the  solution.  The  density  at  t  =  2.3  is  presented  in  Figure  4.9  for 
the  case  and  in  Figure  4.10  for  the  case.  Rectangular  meshes  are  used,  with  four 
different  mesh  sizes  Ax  =  Ay  =  gQ,  40  ’  8o>  respectively. 

We  remark  that  it  is  easy  to  get  negative  density  and/or  pressure  for  this  problem.  In 
both  OTir  P^  and  P^  runs,  we  found  it  necessary  to  perform  a  positivity  correction  procedure. 

For  the  P^  case,  we  check,  for  each  element,  whether  the  density,  as  a  linear  function,  is 
too  close  to  zero  in  the  element.  Speciflcally,  using  the  notation  of  (3.3),  we  check  if 

U  \Ux\  lUj,!  <  “W, 


and,  if  yes,  the  slopes  Ux  and  Uy  are  reduced  by  a  factor 


factor 


1  u 

2  I  Ux  I  "b  I  Uy  I 


The  same  correction  procedure  is  performed  on  the  total  energy.  We  do  not  modify  the  two 
momenta. 

For  the  P^  case,  a  somewhat  stronger  positivity  correction  procedure  is  needed.  We 
check,  for  each  element,  whether  the  density  or  the  total  energy  is  too  close  to  zero.  Using 
the  notation  of  (3.5),  we  check  if 


U  I'Uj,]  lltj,!  \Uxy\  gdUxi]  "h  <  0, 
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Figure  4.4:  Forward  facing  step  problem.  Second  order  P*  results.  Density  p.  30  equally 
spaced  contour  lines  from  p  =  0.090338  to  p  =  6.2365.  Mesh  refinement  study.  Prom  top  to 
bottom:  A.'c  =  Aj/  “  leoj  ^.nd 


24 


Figure  4.5:  Forward  facing  step  problem.  Third  order  results.  Density  p.  30  equally 
spaced  contour  lines  from  p  =  0.090338  to  p  =  6.2365.  Mesh  refinement  study.  From  top  to 
bottom:  Aa;  =  Ar/  =  and 
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Figure  4.6:  Forward  facing  step  problem.  The  triangulations  used  for  results  in  Figures  4.7 
and  4.8.  g  is  the  ratio  between  the  typical  size  of  the  triangles  near  the  corner  and  that 
elsewhere. 
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Figiire  4.7;  Forward  facing  step  problem.  Second  order  P*  results.  Density  p.  30  equally 
spaced  contour  lines  from  p  =  0.090338  to  p  =  6.2365.  THangle  code.  Progressive  refinement 
near  the  comer,  using  the  triangulations  shown  in  Figure  4.6. 


Figure  4.8:  Forward  facing  step  problem.  Third  order  restilts.  Density  p.  30  equally 
spaced  contour  lines  from  p  =  0.090338  to  p  =  6.2365.  Ttiangle  code.  Progressive  refinement 
near  the  corner,  using  the  triangulations  shown  in  Figure  4.6. 


for  either  the  density  or  the  total  energy,  and,  if  yes,  all  the  degrees  of  freedom  except  the 
mean 

'^yy't 

of  all  the  components  (density,  two  momenta,  and  total  energy) ,  are  reduced  by  a  factor, 
which  is  the  smaller  of  the  two  quantities 

u 

“1“  Hh/I  "f*  3(l'^3:a:|  “1“ 

computed  from  the  density  and  the  total  energy. 

We  remark  that  the  positivity  correction  procedures  described  above  are  conservative 
and  do  not  degrade  the  formal  accuracy  of  the  schemes. 

5  Concluding  Remarks 

We  have  presented  the  algorithm  formulation  and  practical  implementation  issues  of  the 
RKDG  (Runge-Kutta  discontinuous  Galerkin)  methods,  for  multidimensional  systems  and 
in  particiilar  for  the  compressible  Euler  equations  of  gas  dynamics.  Numerical  results  are 
shown.  We  conclude  in  particular  that  for  detailed  features  in  the  flow,  such  as  the  structure 
near  the  triple  Mach  stem  in  the  double  Mach  reflection  problem,  a  higher  order  method 
gives  better  cpu  performance  than  a  lower  order  one,  to  obtain  the  same  resolution.  We  also 
conclude  that  triangular  elements  and  rectangular  elements  perform  in  a  similar  way. 
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